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Abstract 


This  report  describes  numerical  predictions  of  hydrodynamic  forces  and  motions  in 
the  time  domain  for  an  unappended  ship  hull  in  waves.  The  hydrodynamic  forces  in 
the  time  domain  are  evaluated  using  previously  computed  values  in  the  frequency 
domain,  which  are  based  on  the  three-dimensional  Green  function  for  a  body  with 
zero  forward  speed.  Forward  speed  effects  are  included  in  both  the  frequency  and 
time  domains  using  suitable  approximations.  Due  to  the  underlying  importance  of 
numerical  predictions  of  hydrodynamic  forces  in  the  frequency  domain,  attention 
has  been  given  to  their  accuracy,  robustness,  and  efficiency.  Comparisons  between 
motion  predictions  in  the  frequency  and  time  domains  show  excellent  agreement, 
suggesting  correct  theoretical  and  numerical  implementations.  Comparisons  of  nu¬ 
merical  heave  and  pitch  predictions  with  model  test  data  for  the  frigate  HALIFAX 
show  very  good  agreement.  Time  domain  simulations  of  ship  motion  execute  sig¬ 
nificantly  faster  than  real  time,  indicating  that  the  approach  will  be  practical  for 
modelling  and  simulation  applications. 


Resume 


Dans  ce  rapport  nous  decrivons  les  predictions  numeriques  calculees  dans  le  do- 
maine  temporel  des  forces  hydrodynamiques  et  des  mouvements  d’une  coque  sans 
appendices  dans  les  vagues.  Nous  avons  evalue  les  forces  hydrodynamiques  dans 
le  domaine  temporel  en  utilisant  les  valeurs  calculees  anterieurement  dans  le  do- 
maine  frequentiel  avec  une  fonction  de  Green  tridimensionnelle  pour  un  corps  de 
vitesse  avant  nulle.  Les  effets  de  la  vitesse  avant  ont  ete  ajoutes  dans  les  domaines 
frequentiel  et  temporel  en  utilisant  des  approximations  adaptees.  Etant  donne  T  im¬ 
portance  intrinseque  des  predictions  numeriques  des  forces  hydrodynamiques  dans 
le  domaine  frequentiel,  nous  avons  apporte  une  attention  particuliere  a  la  precision, 
a  la  robustesse  et  a  l’efficacite  des  calculs.  L’ accord  entre  les  predictions  des  mou¬ 
vements  dans  les  domaines  temporel  et  frequentiel  est  bon,  ce  qui  suggere  que  la 
theorie  et  son  implementation  logicielle  soient  correctes.  L’ accord  entre  les  predic¬ 
tions  numeriques  du  soulevement  et  du  tangage  pour  la  fregate  HALIFAX  avec  les 
donnees  sur  maquette  est  tres  bon.  Le  calcul  des  simulations  des  mouvements  du 
navire  dans  le  domaine  temporel  est  beaucoup  plus  rapide  que  le  calcul  en  temps 
reel,  ce  qui  indique  que  l’on  pourra  utiliser  notre  methode  pour  la  modelisation  et 
les  simulations. 
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Executive  summary 

Introduction 

DRDC  Atlantic  is  developing  a  new  object-oriented  library  for  simulation  of  ship 
motions  in  waves.  A  previous  report  presented  predictions  of  ship  hydrodynamic 
forces  and  motions  in  the  frequency  domain  using  a  panel  method.  This  report 
presents  a  method  for  extending  hydrodynamic  force  and  motion  predictions  to  the 
time  domain  using  previously  computed  values  for  the  frequency  domain. 


Principal  Results 

Improvements  have  been  made  to  accuracy,  robustness,  and  computational  speed 
of  hydrodynamic  force  predictions  in  the  frequency  domain.  New  hydrodynamic 
force  predictions  in  the  time  domain  produce  motions  which  are  in  agreement  with 
predictions  in  the  frequency  domain,  suggesting  correct  theoretical  and  numerical 
implementations.  Comparisons  of  heave  and  pitch  predictions  with  experimental 
results  show  very  good  agreement. 


Significance  of  Results 

New  object-oriented  software  is  now  available  for  simulating  motions  in  waves 
of  an  unappended  ship.  The  new  software  uses  a  panel  method,  ensuring  fidelity 
for  non-slender  hull  forms  such  as  the  Canadian  Navy’s  KINGSTON  class.  The 
inclusion  of  nonlinear  forces  from  buoyancy  and  incident  waves  can  contribute  to 
accuracy  of  ship  motion  predictions  in  higher  sea  states.  Both  linear  and  nonlinear 
predictions  run  significantly  faster  than  real  time,  which  is  important  for  simulation 
applications. 


Future  Plans 

In  the  near  term,  forces  from  rudders,  bilge  keels,  and  other  appendages  will  be 
incorporated  into  numerical  predictions  of  ship  motions  in  waves.  In  the  longer 
term,  the  motions  of  a  freely  maneuvering  ship  in  waves  will  be  predicted. 
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Time  Domain  for  an  Unappended  Ship  Hull;  DRDC  Atlantic  TM 
2003-104;  Defence  R&D  Canada  -  Atlantic. 


DRDC  Atlantic  TM  2003-1 04 


Sommaire 


Introduction 

RDDC  assemble  une  nouvelle  bibliotheque  de  logiciels  orientes  objet  qui  seront  utilises 
pour  simuler  le  mouvement  des  navires  dans  les  vagues.  Dans  un  rapport  precedent,  nous 
avons  presente  les  predictions  des  forces  hydrodynamiques  exercees  sur  le  navire  et  de  ses 
mouvements,  calculees  dans  le  domaine  des  frequences  avec  la  methode  des  panneaux. 
Dans  ce  rapport,  nous  presentons  une  methode  permettant  d’etendre  les  predictions  des 
forces  hydrodynamiques  et  des  mouvements,  au  domaine  temporel,  a  partir  des  valeurs 
deja  calculees  dans  le  domaine  frequentiel. 

Resultats  principaux 

Nous  avons  ameliore  la  precision,  la  robustesse  et  la  vitesse  du  calcul  des  predictions  des 
forces  hydrodynamiques  dans  le  domaine  des  frequences.  Les  nouvelles  predictions  pour  la 
force  dans  le  domaine  temporel  se  traduisent  par  des  mouvements  qui  concordent  avec  les 
predictions  dans  le  domaine  frequentiel,  ce  qui  suggere  que  la  theorie  et  son  implementation 
logicielle  sont  correctes.  L’ accord  entre  les  predictions  pour  le  soulevement  et  le  tangage  et 
les  resultats  experimentaux  est  bon. 

Importance  des  resultats 

Un  nouveau  logiciel  oriente  objet  permettant  de  simuler  les  mouvements  d’un  navire  muni 
d’aucun  appendice  dans  des  vagues  est  maintenant  disponible.  Les  predictions  de  ce  nou¬ 
veau  logiciel  sont  effectuees  a  partir  de  la  methode  des  panneaux,  ce  qui  garantit  l’exacti- 
tude  des  calculs  pour  les  coques  larges  conime  celles  des  navires  de  la  classe  Kingston  de 
la  Marine  canadienne.  L'ajout  des  forces  non  lineaires  de  la  flottabilite  et  des  vagues  inci- 
dentes  pourra  contribuer  a  la  precision  des  predictions  des  mouvements  des  navires  dans  les 
mers  tres  agitees.  Le  calcul  lineaire  et  non  lineaire  des  predictions  est  beaucoup  plus  rapide 
que  le  calcul  en  temps  reel,  un  element  important  pour  les  simulations. 

Travaux  ulterieurs  prevus 

Nous  incorporerons  prochainement  les  forces  exercees  par  les  gouvernails,  les  quilles  de 
roulis  et  autres  appendices  aux  predictions  numeriques  des  mouvements  des  navires  dans 
les  vagues.  A  plus  long  terme,  on  pourra  predire  des  mouvements  de  navires  manceuvrant 
librement  dans  les  vagues. 

Kevin  McTaggart;  2003;  Hydrodynamic  Forces  and  Motions  in  the 
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1  Introduction 


Simulations  of  ship  motions  in  waves  are  required  for  numerous  naval  applications. 
When  evaluating  performance  of  crew  and  ship  systems,  the  influence  of  ship  mo¬ 
tions  is  often  significant.  The  safety  of  a  ship  will  be  highly  dependent  on  its  ability 
to  resist  capsize  and  wave-induced  structural  loads. 

To  facilitate  simulation  efforts,  DRDC  Atlantic  is  developing  a  new  object-oriented 
library  for  predicting  ship  motions  in  waves.  During  development  of  this  library, 
attention  is  being  given  to  balancing  requirements  for  fidelity,  efficiency,  and  ro¬ 
bustness.  It  was  decided  that  these  requirements  could  be  satisfied  by  computing 
ship  motions  in  the  time  domain  using  hull  hydrodynamic  forces  derived  from  com¬ 
putations  in  the  frequency  domain.  Reference  1  describes  computation  of  hull  hy¬ 
drodynamic  forces  in  the  frequency  domain  using  a  panel  method.  The  present 
report  describes  the  evaluation  of  hydrodynamic  forces  in  the  time  domain  using 
previously  computed  results  in  the  frequency  domain.  Section  2  gives  a  review  of 
hydrodynamic  coefficients  in  the  frequency  domain  based  on  Reference  1.  The  the¬ 
ory  for  obtaining  time  domain  forces  from  frequency  domain  values  is  presented 
in  Section  3.  Section  4  gives  particulars  for  the  ship  HALIFAX,  which  is  used  for 
sample  computations  throughout  the  remainder  of  the  report.  Improvements  to  fre¬ 
quency  domain  computations  are  given  in  Section  5,  followed  by  sample  frequency 
domain  coefficients  in  Section  6.  Section  7  gives  sample  time  domain  coefficients 
based  upon  previously  computed  frequency  domain  values.  Section  8  discusses 
computation  of  ship  motions  in  the  time  domain  using  time  domain  hydrodynamic 
coefficients.  To  verify  the  new  time  domain  approach.  Section  9  compares  mo¬ 
tions  computed  in  the  time  and  frequency  domains.  The  extension  of  time  domain 
predictions  to  include  nonlinear  contributions  from  buoyancy  and  incident  waves  is 
described  in  Section  10,  followed  by  validation  of  computations  using  experimental 
data  in  Section  11.  Section  12  gives  concluding  remarks. 
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2  Review  of  Hydrodynamic  Coefficients  in 
the  Frequency  Domain 


As  mentioned  previously,  this  report  describes  prediction  of  hull  hydrodynamic 
forces  in  the  time  domain  using  previously  computed  frequency  domain  values. 
For  completeness,  a  review  of  frequency  domain  coefficients  from  Reference  1  is 
presented  here. 

Figure  1  shows  the  coordinate  system  used  for  ship  motion  computations.  The  ship 
is  assumed  to  move  at  steady  forward  speed  and  heading,  and  the  coordinate  system 
is  a  translating  earth  system  moving  with  the  steady  forward  speed.  The  origin  of 
the  translating  earth  system  is  located  at  the  ship  centre  of  gravity  when  the  ship  is 
in  calm  water.  The  location  of  the  translating  earth  system  is  not  affected  by  wave- 
induced  ship  motions.  Note  that  the  present  coordinate  system  is  slightly  different 
from  that  used  for  solution  of  hydrodynamic  forces  in  Reference  1,  for  which  the 
vertical  z  axis  passed  through  the  centre  of  gravity  but  had  its  origin  at  the  calm 
water  surface.  Reference  1  used  a  different  origin  because  many  computational 
terms  were  dependent  upon  elevation  relative  to  the  calm  water  surface.  In  the 
present  axis  system,  elevation  relative  to  the  calm  water  surface  is  given  by: 

zwt  =  z  +  (1) 

where  is  the  elevation  of  the  ship  centre  of  gravity  relative  to  the  calm  water 
surface  when  the  ship  is  at  rest.  Figure  2  gives  the  sea  direction  convention,  with  a 
heading  (3  of  180  degrees  representing  head  seas. 

Oscillatory  ship  motions  in  the  time  domain  are  related  to  ship  motions  in  the  fre¬ 
quency  domain  as  follows: 

r)j(t)  =  Really  exp(i  uie  t)}  for  j  =  1  —  6  (2) 

where  cue  is  wave  encounter  frequency  given  by: 

c oe  =  \uj  —  ki  U  cos/3|  (3) 

where  tui  is  incident  wave  frequency,  ki  is  incident  wavenumber,  U  is  forward  ship 
speed,  and  f3  is  incident  wave  direction.  Water  depth  is  assumed  to  be  large  (i.e., 
greater  than  half  the  incident  wavelength),  and  the  following  dispersion  relation 
applies: 


ki  = 


9 


(4) 
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Figure  1:  Coordinate  System  for  Solution  of  Hydrodynamic  Forces 


Figure  2:  Sea  Direction  for  Solution  of  Hydrodynamic  Forces 
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where  g  is  gravitational  acceleration.  The  equations  of  ship  motions  in  the  fre¬ 
quency  domain  are  given  by: 


{-^([M]  +  [A])  +  iu.[B]  +[C]}  in]  =  {F‘  +  Fd}  (5) 


where  [M]  is  the  ship  mass  matrix,  [A]  is  the  added  mass  matrix,  [B]  is  the  damping 
matrix,  [C\  is  the  stiffness  matrix,  {  F1 }  is  the  incident  wave  excitation  vector,  and 
{Fd}  is  the  diffracted  wave  excitation  vector.  The  added  mass  and  damping  matrix 
terms  are  given  by: 


Ajk 

Bjk 


P_ 

Lde 


*sb  L 


Imag  ((/>fc)  + 


U 

—  Real 

cce 


-P 


>sb  L 


Real(</>* 


—  Imag 

(jJe 


rij  dS  (6) 

rij  dS  (7) 


where  p  is  water  density,  Sb  denotes  the  wetted  surface  of  the  ship,  cpk  is  the  complex 
velocity  potential  due  to  oscillatory  motions  for  mode  k,  and  rij  is  the  outward 
normal  component  for  mode  j.  Note  that  Equation  (6)  has  a  correction  to  a  sign 
error  for  the  speed  term  in  Reference  1 .  The  above  equations  can  be  rewritten  to 
show  speed  dependent  terms  as  follows: 

Ajk  =  —  —  f  lmag(0fc)ni  dS  -  —  —  [  Real  rij  dS  (8) 

Jsb  Ue  UJe  JSb  V  dx  J 

Bjk  =  -p  J  Real(0fc)  n3  dS  +  P  ~  J  Ima§  n'J  dS 


At  zero  and  infinite  frequency  limits,  the  real  part  of  the  velocity  potential  goes 
to  zero,  and  the  velocity  potential  can  be  represented  in  terms  of  normalized  real 
quantities  as  follows: 


0( ue)  =  i  4>\0)  uje  foro;e  — >  0  (10) 

(f>( ue)  =  i  <f>'(oo)  uje  for  uje  — >  oo  (11) 

At  these  limits  the  added  mass  and  damping  can  be  evaluated  as: 


Ajk 

=  ~P 

f  fik'nj  dS 

(12) 

>sb 

Bjk 

=  pU 

its* 

(13) 

For  pitch  and  yaw,  the  velocity  potentials  at  forward  speed  include  contributions 
from  heave  and  sway  respectively  as  follows: 

( h(U,ue )  =  05(O  ,ue)  +  —  03(O,cce)  (14) 

X  UJe 

( j>e(U,LUe )  =  06(O  ,Ue)  -  -  02(O,CUe)  (15) 

l  L Of, 
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Added  mass  and  damping  coefficients  for  pitch  can  be  expressed  as: 


A/5  — 


P 


P  U 


- /  Imag  (<^5(0,  u>e))  rij  dS  + - /  Real  (cf>3)njdS 

^ e  J  Sb  ^ e  J  Sb 

/  Real  (d^0^ 


^ e  W  e  J  Sb 
p  (V 


Ue  \  U, 


V  dx 

Imag 


rij  dS 


>sb 


{  ^tA 

\dx  J 


rij  dS 


U 


(16) 


Bj5  =  —p  /  Real(05(O,  ue))  rij  dS  —  p —  /  Imag  (03)  rij  dS 


Jsb 

P  —  I  Imag 


uj, 


u. 


~  P 


e  JSb 

u_ 

c oP 


(905  (0  ,ue 


Real 


'sb 


dx 
<903  \ 


e  JSb 


•rij  dS 


(17) 


At  the  zero  and  upper  frequency  limits,  the  above  equations  can  be  re-written  as: 

2  /'  drhi 

n,  dS  (18) 


Aj5  =  -p  ( 05 )  rij  dS  -  p  (Jpj 


Bj5  =  -pU  (0£) rij  dS  +  pU 


'sb 


'sb 


Jsb 

M. 

dx  Uj 


dx  '  ] 
n,  dS 


(19) 


For  yaw  motion,  added  mass  and  damping  terms  are: 


A?  6  — 


P 


P  U 


- /  Imag  (06 (0,  uje))  rij  dS  — - /  Real  (02)  rij  dS 

J  Sb  ^e.  J  Sb 


-  E  f  V  » 

UJe  ue  Jsb  V  dx 


( dMhrA]  n,  dS 


P_  (U_ 

^ e  \^e 


Imag 


'sb 


(  AtA 
\dx  J 


Bj6 


■rij  dS 


U 


(20) 


=  —p  /  Real(06(O,u/e))  rij  dS  +  p —  /  Imag  (02)  rij  dS 


Jsb 

P  —  I  Imag 


O’, 


P 


JSb 
f/X  2 

C Op 


006 (0,  U/e 


Real 


'Sb 


dx 

<902  , 
a7 )  % 


e  JSb 


rij  dS 


Yaw  motion  terms  at  the  zero  and  upper  frequency  limits  are: 

2 


Aj 6  =  -P  J  (06 )nj  dS  +  p 

BjQ  =  pU  [  (02 )nj  dS  +  p  U 


9  A  Ao 

— — n?  ciO 
ax 


'Sb 


Jsb 

<906 

— —  n?-  do 
ax 


(21) 

(22) 

(23) 
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For  computational  purposes,  it  is  useful  to  be  able  to  express  ship  hydrodynamic 
coefficients  at  forward  speed  as  functions  of  zero  speed  coefficients  as  follows: 


Ajk(U,  ^Je) 

A-jk (0}  ^Je) 

U  dBjk(0,uje) 
l dx 

for  k  =  1  —  4 

(24) 

Bjk(U,  e ) 

-^jk  (0?  ^e) 

^  dAjk(0,ue) 
dx 

for  k  =  1  —  4 

(25) 

AjS {JJ i  ^e) 

(0)  ^e) 

U 

—  Bj3(0,ue)  + 

U  <9-Bj5(0,  ue) 

ul  dx 

*  ©‘ 

dAj3{Q,ue) 

(26) 

dx 

Bj5{U,ue)  =  Bj5(0,cue)  +  U  Aj3(0,u>e)  -  U  dA^u^ 
fU_V  dBj3(0,uje) 

oe)  dx 

Aje(U,Lue)  =  AjQ(0,Lue )  +  —  Bj2(0,u>e)  +  — 


c ot 

U\2  dAj2(0,me) 

LUeJ  dx 

Bje(U,LUe)  =  Bj6(0,ue)  —  U  Aj2(0,coe)  —  U 

U\2  <9-5)2  (0,  X>e) 

,  LUe  J  dx 


dx 


dAj6(0,ue 

dx 


(27) 


(28) 


(29) 


The  above  equations  are  the  basis  for  a  ship  hydrodynamic  coefficient  database 
which  gives  coefficients  at  arbitrary  forward  speed  based  on  stored  coefficients  for 
zero  forward  speed. 


3  Hydrodynamic  Forces  for  Transient 
Motions 


The  hydrodynamic  forces  for  transient  motions  of  floating  bodies  are  discussed  in 
several  references.  Mei  [2]  and  Wehausen  [3]  give  comprehensive  overviews.  Mc- 
Taggart  [4]  discusses  transient  motions  of  icebergs  during  drifting  and  collisions 
with  offshore  structures.  Magee  [5]  gives  a  good  overview  of  transient  hydrody¬ 
namic  forces  acting  on  ships,  and  includes  forward  speed  effects.  Fonseca  and 
Guedes  Soares  [6]  present  theory  and  results  for  time  domain  ship  motions  based 
on  frequency  domain  coefficients. 

For  a  floating  body  undergoing  transient  motions,  the  equations  of  motion  can  be 
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expressed  as: 


m  +  im  oo)])  w))  +  im\  {<?w} 

+  [  [K{U,t- t)]{7)(t)}  dr 

J  — OO 

+  ([Cl  +  [c(t/)[)  {»,(«)}  =  {F'(t)  +  FD(t)}  (30) 

where  [b(U)]  is  the  speed  dependent  damping  matrix  in  the  time  domain,  [K (U.  t  —  t 
is  the  retardation  function  matrix  for  delay  time  t  —  r,  and  [c(U)]  is  the  speed  de¬ 
pendent  stiffness  matrix  in  the  time  domain.  Terms  in  the  above  equation  can  be 
evaluated  using  results  from  frequency  domain  computations. 


3.1  Retardation  Functions 

The  retardation  functions  Kjk(U,T )  can  be  determined  from  frequency  domain 
added  mass  or  damping  coefficients  as  follows: 

2  f°° 

Kjk{r)  = - /  [Ajk(U,Lue)  -  Ajk(U,  oo)]  ue  sino;eT  due  (31) 

71  Jo 

2 

Kjk{r)  =  -  Bjk(U,uje)  costUeT  due  (32) 

77  Jo 

When  selecting  between  the  above  two  equations  for  evaluating  retardation  func¬ 
tions,  it  should  be  noted  that  damping  coefficients  approach  their  infinite  frequency 
limits  much  more  quickly  than  do  added  masses.  Given  the  problems  with  irregular 
frequencies  that  can  occur  at  higher  frequencies,  it  was  decided  to  compute  re¬ 
tardation  functions  based  on  damping  where  possible.  Fonseca  and  Guedes  Soares 
present  a  similar  approach  based  on  hydrodynamic  coefficients  evaluated  using  strip 
theory. 

It  is  convenient  to  express  retardation  functions  in  terms  of  speed-independent  and 
speed-dependent  terms.  For  modes  for  which  velocity  potentials  are  assumed  to 
be  independent  of  ship  speed  (i.e.,  surge,  sway,  heave,  and  roll),  the  retardation 
functions  can  be  expressed  as: 


Kjk(U,  t)  =  K°jk(r)  +  UK%{t)  for  k  =  1-4  (33) 

The  speed-independent  retardation  functions  are  evaluated  using: 

2  f°° 

Jv°  (r)  —  —  Bjk(0,LUe )  cosu;eT  doje  for  k  =  1  —  6  (34) 

71  Jo 
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Note  that  the  above  equation  is  valid  for  all  six  degrees  of  freedom.  The  speed 
dependent  terms  based  on  Equations  (25)  and  (32)  are: 


Kfk(r)  = 


n 


poo 

~dAjk(0,Loe) 

dAjk(0,oo)~ 

L 

dx 

dx 

cosc Oft  duP 


for  k  —  1  —  4 


(35) 


For  pitch  and  yaw,  the  velocity  potentials  are  influenced  by  forward  speed,  and  the 
retardation  functions  are  expressed  as: 

Kjk(U,r)  =  K°k{r)  +  U  R%(r)  +  U2  A^r)  for  A;  =  5  —  6  (36) 

The  speed-dependent  coefficients  due  to  pitch  motion  are  evaluated  by: 


2 
7 r 
2 
71 


dAj5(0,ue)  dAj  5(0,oo) 


dx 


dx 


COS  LUeT  duje 


UP 


2 

7T 


dBjs(0,  tue)  cos u>eT 


dx 


duje 

(37) 

du)e 

(38) 

When  evaluating  in  Equation  (37),  the  contribution  from  the  heave  velocity  po¬ 
tential  is  evaluated  using  heave  damping  rather  than  added  mass  because  it  has  been 
found  to  give  better  numerical  accuracy.  Both  Equations  (37)  and  (38)  have  inte¬ 
grands  with  encounter  frequency  c oe  or  its  square  in  the  denominator.  Fortunately, 
these  integrands  are  not  problematic  as  cue  approaches  zero  because  damping  is  zero 
at  zero  encounter  frequency. 


Similarly,  the  speed-dependent  coefficients  due  to  yaw  motion  are  evaluated  by: 


2 

71 

2 

7T 

2 

71 


'0 


dAj6(0,u;e)  cL4j6(0,  00) 


dx 


dx 


COS  U)pT  dUp 


/  D  fn  S  sm  ueT 

/  Bj2{0,uje)  - due 

1 0 
roo 


dBj2(0,LUe)  COS  LUeT 


dx 


dujp 


out 


(39) 

(40) 


3.2  Hydrodynamic  Damping  and  Stiffness 

Hydrodynamic  damping  and  stiffness  terms  are  related  to  frequency  domain  values 
as  follows: 


bjk(U) 


jj  dAjk(0}  00) 
dx 


for  k  =  1  —  4 


(41) 
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(42) 


WEO 

Cjs(^) 

Cj6(£/) 


_u»^06)  +  p^(0jCo) 


=  -  [/ 


=  -u[ 


dx 

dAj6(Q,  oo) 
dx 

2  <9^3(0,  0) 


-  U  Aj2( 0,  00) 


=  u 


dx 

2  <9^2(0, 0) 
dx 


(43) 

(44) 

(45) 


Note  that  the  two  stiffness  terms  arise  from  converting  frequency  domain  added 
mass  terms  to  stiffness  terms  using  multiplication  by  —  u2e.  The  stiffness  terms  are 
based  on  sway  and  heave  added  masses  at  zero  frequency. 


3.3  Wave  Excitation  Forces 

Linear  wave  excitation  forces  in  the  time  domain  can  be  evaluated  using  the  fre¬ 
quency  domain  forces  presented  in  Reference  1.  For  a  linear  wave  system,  the 
incident  wave  elevation  in  the  time  domain  can  be  expressed  as: 

( i(x,y,t )  =  a  cos  [uet  +  ej  —  ki  (x  cos  (3  —  ysin/3)]  (46) 

where  a  is  incident  wave  amplitude,  ej  is  the  phase  of  the  incident  wave  crest  (posi¬ 
tive  for  time  lead)  at  the  origin,  kj  is  the  incident  wavenumber,  and  / 3  is  the  incident 
wave  direction.  Wave  excitation  forces  in  the  time  domain  can  be  expressed  in 
terms  of  frequency  domain  values  as  follows: 

Fj  (t)  =  |F7(f/,/3,^)|  cosset  +6/  +  e  (F1  (U,  (3,^))]  (47) 

FP(t)  =  \Fd(U,P,U!)\  cos[uet  +6!  +  e(FD(U,P,ui))]  (48) 

where  F1  (U,  /3,ui)  is  the  complex  incident  excitation  force  in  the  frequency  do¬ 
main  and  F,)(U.  /3,ui)  is  the  complex  diffraction  excitation  force  in  the  frequency 
domain.  The  phases  of  the  incident  and  diffracted  excitation  forces  for  incident 
waves  with  zero  phase  angles  are  given  by: 


e(  F<(U,f3,uj,)) 
e(FD(U,/3,  a.,)) 


arctan 

arctan 


Imag  {F1  (U,  f3,u>i)} 
Real  {FPU^,^)} 
Imag  {Fd( C7, /3,a;/)} 
Real  {FD(U,P,uT)} 


(49) 

(50) 
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4  CPF  Hydroelastic  Model  Hull  for  Sample 
Computations  of  Hydrodynamic 
Coefficients 


To  demonstrate  results  of  computations  this  report  uses  the  HALIFAX  class,  also 
known  as  the  Canadian  Patrol  Frigate  (CPF).  To  enable  comparison  of  motion  pre¬ 
dictions  with  model  test  results,  the  loading  condition  for  the  present  computations 
is  the  deep  departure  condition  for  the  CPF  hydroelastic  model,  as  described  in 
Reference  7.  Note  that  a  slightly  different  loading  condition  was  used  for  sample 
computations  of  frequency  domain  coefficients  in  Reference  1 . 

Table  1  gives  the  main  particulars.  Three  different  hydrodynamic  meshes  have  been 
generated,  as  shown  in  Figures  3  to  5.  The  panel  colours  in  each  figure  represent 
depth  below  the  waterline.  Table  2  gives  hydrostatic  properties  computed  using 
each  mesh.  Preliminary  hydrodynamic  computations  indicated  that  the  x  deriva¬ 
tive  of  added  mass  was  somewhat  sensitive  to  the  limiting  aspect  ratio  (lenght- 
wise/girthwise)  used  for  panelling  the  hull.  Limiting  aspect  ratios  of  3  and  5  were 
tested,  with  a  value  of  3  being  ultimately  chosen  due  to  better  convergence  of  results 
with  increasing  number  of  panels. 

Table  1:  Main  Particulars  for  HALIFAX  Class  Frigate,  CPF  Hydroelastic  Model 
Deep  Departure  Condition 


Length,  L 
Beam,  B 

Midships  draft,  Tmid 
Trim  by  stern,  ts 
Displacement,  A 
Vertical  centre  of  gravity,  KG 
Dry  roll  radius  of  gyration  rxx 
Dry  pitch  radius  of  gyration  ryy 
Dry  yaw  radius  of  gyration  rzz 


124.5  m 

14.8  m 
4.97  m 
-0.04  m 

4655  tonnes  (fresh  water) 
6.26  m 
5.82  m 

28.8  m 
28.8  m 
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Figure  3:  Coarse  Mesh  for  CPF  Hydroelastic  Model 


Figure  4:  Medium  Mesh  for  CPF  Hydroelastic  Model 
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Figure  5:  Fine  Mesh  for  CPF  Hydroelastic  Model 


Table  2:  Hydrostatic  Properties  for  Different  Panel  Meshes 


Coarse 

Medium 

Fine 

Nominal  panel  size  (m2) 

5.0 

2.5 

1.0 

Number  of  panels  on  port  side  N^ort 

224 

433 

1040 

Volume  (m3) 

4465.8 

4488.8 

4500.6 

LCB  aft  of  FP  (m) 

64.663 

64.665 

64.678 

CB  wrt  waterline  (m) 

-1.882 

-1.889 

-1.893 

Wetted  surface  area  (m2) 

1965.6 

1969.4 

1971.2 
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5  Improvements  to  Computations  of 
Frequency  Domain  Coefficients 


The  present  approach  for  computing  time  domain  hydrodynamic  coefficients  re¬ 
quires  accurate  and  efficient  computation  of  frequency  domain  coefficients.  This 
section  presents  several  improvements  that  have  been  made  to  computations  of  fre¬ 
quency  domain  coefficients  subsequent  to  the  work  reported  in  Reference  1. 

5.1  Maintenance  of  Points  on  Waterline  During 
Panel  Transformation 

The  wetted  surface  of  a  hull  is  represented  using  triangular  and  quadrilateral  panels. 
As  discussed  in  Reference  1,  it  is  usually  necessary  to  make  slight  modifications  to 
points  of  quadrilateral  panels  to  ensure  that  all  points  lie  in  the  same  plane.  Pan¬ 
elling  computations  have  been  changed  such  that  a  panel  point  located  on  the  hull 
waterline  remains  on  the  waterline  zwi  =  0,  with  only  x  and  y  coordinates  being 
changed  to  obtain  a  planar  panel.  This  modification  prevents  numerical  inconsis¬ 
tencies  which  can  occur  when  panel  points  are  located  above  the  waterline. 

5.2  Sign  Correction  to  Speed  Dependent  Added 
Mass 

As  stated  in  Section  2,  Equation  9  of  Reference  1  contained  a  sign  error  for  the 
speed  dependent  term  of  added  mass.  Equation  (6)  gives  the  correct  added  mass. 


5.3  Improved  Accuracy  for  x  Derivatives  of  Velocity 
Potential 


Improved  accuracy  has  been  introduced  for  evaluating  the  x  derivatives  of  velocity 
potential.  The  following  equation  is  still  used  for  computing  the  x  derivatives: 


©  ■  (»i  <■> 


(51) 


where  [H]  is  the  influence  matrix  and  {cr}  is  the  vector  of  source  strengths.  The 
influence  matrix  coefficients  were  previously  evaluated  as  follows: 


Hjk  &jk  2  A/: —j  T  (1 


5jk)  An 


•sk 


(9(1  /R{xj,xs)) 
dx? 


dS 
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Hjk 


+  J f  d(l/Ri(xj,xs)) 

4n  JSk  ~  DxSj 

+  -  f  aSf^  dS 

4k  Jsk  oxSj 

1 

^ jk  2  ^x—j  4~  (1  5 jk ) 

_ 1_  /'  (9(l/i?i(fj,fs)) 

4tt  Jsk  dxSj 

+  J_  f  dGoojx^Xs)  dS 

4k  Jsk  dxSj 


dS 

for  CUe  <  w* 

r  (9(1  /R(xj,xs)) 

4k  Jsk  dx3j 

dS 

for  LUe  >  u>l 


dS 


(52) 


(53) 


where  nr_  ?  is  the  x  component  of  the  outward  normal  on  panel  j,  5jk  is  the  Kroe- 
necker  delta  function,  Sk  is  the  surface  of  source  panel  k,  R(xj,xs )  is  the  radius 
from  field  point  Xj  to  source  panel  point  xs,  R\  is  the  radius  from  the  field  point 
to  the  image  of  the  source  panel  point,  Go  is  the  frequency  dependent  Green  func¬ 
tion  relative  to  the  zero  frequency  Green  function,  a is  the  transition  encounter 
frequency  for  indicating  which  form  of  the  frequency  dependent  Green  function 
to  use,  and  G ^  is  the  frequency  dependent  Green  function  relative  to  the  infinite 
frequency  Green  function.  The  selection  of  which  of  the  above  equations  to  use 
depends  upon  a  specified  transition  encounter  frequency  a For  the  case  j  =  k 
considering  the  flow  from  a  source  panel  on  itself,  the  above  equations  are  based 
on  the  following  assumption: 


'sk 


d(l/R(xj,xs)) 

6Xr,: 


dS 


-2  7r  nx-j 


for  j  =  k 


(54) 


Although  this  equation  is  a  good  approximation,  it  has  been  found  that  it  should  not 
be  used  because  computed  x  derivatives  of  velocity  potentials  can  be  very  sensi¬ 
tive  to  small  differences  in  influence  coefficients  Hjk.  Consequently,  the  influence 
coefficients  are  now  evaluated  by: 


Hjk 


Hjk 


J_  r  d(l/R(xj,xa)) 

4k  JSk  dxSj 

+  J _  f  dG0(xj,xa ) 
4k  JSk  dxSj 
J_  k  d(l/R(xj,xs)) 
4k  JSk  dxSj 

+  J_  f  dGoo(xj,  xs) 

4k  JSk  dxSj 


1  f  djl/R^Xj.Xs)) 
+  4v r  4  dxSj 

dS  for  lu£  <  u>l 

dS _ 1  f  djl/R^XjtXs)) 

4k  JSk  dxSj 

dS  for  LUe  >  u>l 


dS 


dS 


(55) 


(56) 
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5.4  Gauss  Quadrature  for  Integration  from  Source 
Panel 

When  integrating  the  Green  function  from  a  source  panel,  the  frequency  dependent 
portion  was  previously  approximated  using  the  centroid  of  the  source  panel.  A  new 
option  has  been  introduced  for  using  Gauss  quadrature  (see  Cook  [8])  on  source 
panels  for  integrating  the  frequency  dependent  portion  of  the  Green  function.  The 
new  implementation  has  options  for  using  1,  4,  or  9  quadrature  points  per  source 
panel. 

5.5  Galerkin  Method  for  Evaluation  of  Velocity 
Potential  Influence  Coefficients 

Previous  computations  of  influence  matrices  were  based  on  the  velocity  potential  at 
the  centroid  of  each  panel.  The  Galerkin  approach  allows  influence  matrix  terms  to 
be  based  on  weighted  values  of  the  Green  function  evaluated  at  1,  4,  or  9  quadrature 
points  on  each  field  panel.  Sclavounos  and  Lee  [9]  give  a  detailed  discussion  of  this 
approach. 

The  accuracy  of  computed  velocity  potentials  will  increase  with  the  number  of 
quadrature  points  used  with  the  Galerkin  method  and  for  integrating  the  frequency 
dependent  portion  of  the  Green  function  from  source  panels.  This  increased  accu¬ 
racy  comes  at  the  expense  of  higher  computation  times  due  to  larger  numbers  of 
computations  of  Green  functions.  Figure  6  shows  heave  damping  and  added  mass 
x  derivative  for  two  different  sets  of  computations.  The  first  set  of  computations 
uses  4  quadrature  points  per  panel  applied  to  the  Galerkin  method  and  integration 
of  the  frequency  dependent  term  from  source  panels.  The  second  set  of  computa¬ 
tions  is  based  on  centroid  approximations.  At  higher  frequencies  (the  first  irregular 
frequency  and  higher),  the  numerical  results  are  sensitive  to  accuracy  of  influence 
coefficients;  thus,  the  higher  order  computations  give  better  results,  particularly  for 
the  x  derivative  of  added  mass. 


5.6  Additional  C  Code  for  Faster  Computations 

As  discussed  in  Reference  1,  most  of  the  present  numerical  predictions  have  been 
developed  in  Python  [10],  a  high-level,  interpreted  computer  language.  When  de¬ 
veloping  Python  applications,  it  is  common  practice  to  determine  which  portions 
are  CPU  intensive  and  to  re-write  them  in  C  to  obtain  the  performance  of  compiled, 
optimized  code.  Previously,  the  frequency  dependent  portion  of  the  Green  func¬ 
tion  was  written  in  C.  New  C  code  has  been  developed  to  evaluate  the  frequency 
independent  terms  of  the  Green  function.  In  addition,  a  new  C  function  has  been 
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Figure  6:  Dimensionless  Forces  from  Zero  Speed  Heave  Damping  and  x 
Derivative  for  HALIFAX  at  20  Knots,  Medium  Mesh,  with  Centroid  and  Four  Point 
Approximations 
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developed  to  evaluate  influence  matrices  for  the  frequency  dependent  portion  of  the 
Green  function  given  input  vectors  of  panel  points. 


5.7  Removal  of  Irregular  Frequencies  from  Database 

Computations  of  hydrodynamic  coefficients  can  be  adversely  influenced  by  irreg¬ 
ular  frequencies,  at  which  mathematical  solutions  of  velocity  potentials  are  incor¬ 
rect.  Me  Taggart  [11]  and  Ando  [12]  discuss  the  occurrence  of  irregular  frequencies 
in  detail. 


For  modelling  and  simulation  work,  it  is  likely  that  hydrodynamic  coefficients  will 
typically  be  determined  using  interpolation  from  a  previously  generated  database 
of  hydrodynamic  coefficients.  The  accuracy  of  interpolated  coefficients  will  be  im¬ 
proved  if  the  database  does  not  have  hydrodynamic  coefficients  adversely  affected 
by  irregular  frequencies.  Fortunately,  a  method  is  available  for  indicating  likely  ir¬ 
regular  frequencies  which  should  be  excluded  from  the  database.  As  was  presented 
in  Reference  1,  radiation  source  strengths  on  the  port  side  of  a  symmetrical  hull  can 
be  solved  by  satisfying  the  following: 

[Deven]  {aT-ort}  =  {d^/1 rt/dn}  for  j  =  1,  3,  5  (57) 

[Dodd]  {afrt}  =  {d^/dn}  for  j  =  2,4,6  (58) 

where  [Deven]  is  the  influence  matrix  for  normal  velocity  on  the  hull  port  side  due 

to  sources  on  the  port  side  when  velocity  potential  is  an  even  function  of  y,  } 

is  the  vector  of  source  strengths  on  the  port  side  for  motion  mode  j,  {<9 (j)p°rt /dn} 
is  the  vector  of  hull  normal  velocities  on  the  port  side,  and  is  the  influence 

matrix  for  normal  velocity  on  the  hull  port  side  due  to  sources  on  the  port  side  when 
velocity  potential  is  an  odd  function  of  y.  The  normal  velocity  influence  matrices 
are  given  by: 


^j-yeven' 


^j-yport,port  j  _|_  ^j-yport,star  j 
^j-yport,port  j  _  ^^portjStarj 


(59) 

(60) 


where  [£)P°ri>Port]  is  the  influence  matrix  for  normal  velocities  on  the  port  side  due  to 
sources  on  the  port  side  and  is  the  influence  matrix  for  normal  velocities 

on  the  port  side  due  to  sources  on  the  starboard  side. 


As  suggested  by  Lee  and  Sclavounos  [13],  irregular  frequencies  are  indicated  by 
high  condition  numbers  for  the  influence  matrices  |  D,  ven  and  [Z}ocM]  for  their  re¬ 
spective  motion  modes.  The  condition  number  for  a  matrix  [D]  can  be  computed  as 
follows  [14]: 


a  =  IMI  X  || [zr1] 


(61) 
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where  1 1  [D]  \  |  is  the  norm  of  matrix  [D]  and  1 1  [D]  1 1 1  is  the  norm  of  its  inverse.  The 
norm  of  [D\  is  given  by: 


N 

1 1 [D] 1 1  =  max y;  \Djk\  (62) 

k= 1 

and  the  norm  of  [D]^1  is  evaluated  in  a  similar  manner.  The  most  computation¬ 
ally  intensive  aspect  of  using  the  above  equation  is  evaluation  of  the  inverse  [D-1]. 
Fortunately,  this  is  not  a  big  obstacle  when  evaluating  hydrodynamic  coefficients 
because  the  number  of  degrees  of  freedom  is  typically  less  than  2000.  Further¬ 
more,  the  computed  inverse  can  be  used  for  computing  source  strengths  using  direct 
matrix  multiplication. 

Figures  7  to  12  show  normalized  hydrodynamic  coefficients  and  matrix  condition 
numbers  for  the  HALIFAX  using  coarse,  medium,  and  fine  panel  meshes.  The 
figures  demonstrate  that  local  maxima  of  matrix  condition  numbers  can  be  indica¬ 
tors  of  irregular  frequencies.  To  reduce  irregular  frequency  effects,  envelopes  of 
maximum  allowable  matrix  condition  numbers  can  be  provided  as  input  for  numer¬ 
ical  computations.  For  a  given  hull  panel  mesh,  separate  envelopes  are  provided 
for  longitudinal  and  lateral  modes.  If  the  computed  matrix  condition  number  for 
a  given  encounter  frequency  exceeds  the  specified  envelope,  then  computations  for 
that  encounter  frequency  are  not  used.  To  determine  appropriate  condition  number 
envelopes,  initial  computations  of  hydrodynamic  coefficients  and  condition  num¬ 
bers  are  required,  followed  by  examination  of  results  such  as  those  given  in  Figures 
7  to  12.  Note  that  the  condition  number  envelopes  will  depend  upon  the  underwater 
hull  geometry,  and  to  a  lesser  extent  upon  the  hull  panel  size  used. 
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Figure  7:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Coarse  Mesh,  Longitudinal  Modes 
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Figure  8:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Coarse  Mesh,  Lateral  Modes 
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Figure  9:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Medium  Mesh,  Longitudinal  Modes 
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Figure  10:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Medium  Mesh,  Lateral  Modes 
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Figure  1 1:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Fine  Mesh,  Longitudinal  Modes 
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Figure  12:  Normalized  Added  Mass  and  Damping,  and  Normal  Velocity  Condition 
Number  for  HALIFAX  Fine  Mesh,  Lateral  Modes 
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5.8  Direct  Evaluation  of  Incident  Wave  Forces  in 
Database  of  Radiation  and  Wave  Excitation 
Forces 

Reference  1  introduced  a  database  class  of  radiation  and  wave  excitation  forces 
in  the  frequency  domain.  This  database  stored  computed  values  for  added  mass, 
damping,  incident  wave  force,  and  diffracted  wave  force.  Using  the  stored  values, 
added  mass  and  damping  could  then  be  quickly  evaluated  through  linear  interpola¬ 
tion  based  on  encounter  frequency.  Similarly,  incident  and  diffracted  wave  forces  in 
unit  amplitude  waves  could  be  quickly  evaluated  through  linear  interpolation  based 
on  ship  speed,  relative  sea  direction,  and  wave  frequency. 

The  database  class  has  been  revised  such  that  incident  wave  forces  are  now  eval¬ 
uated  directly  rather  than  obtained  by  interpolation  of  stored  forces.  Direct  evalu¬ 
ation  of  incident  wave  forces  requires  relatively  little  computation  time;  thus,  this 
revision  increases  accuracy  at  very  little  computational  expense.  Furthermore,  in¬ 
cident  wave  forces  are  usually  greater  in  magnitude  than  diffracted  wave  forces,  so 
the  improvement  provides  a  significant  increase  in  the  accuracy  of  wave  excitation 
forces. 


5.9  Improved  Efficiency  for  Evaluating  Diffracted 
Wave  Potentials 


A  new  computational  approach  has  been  introduced  for  evaluating  diffracted  wave 
potentials.  Recall  that  diffraction  source  strengths  potentials  are  solved  as  follows: 

<">>  =  -  {tr}  (63) 

where  [D]  is  the  influence  coefficient  matrix  relating  hull  normal  velocity  to  source 
strength,  { a D }  is  the  vector  of  hull  diffraction  source  strengths,  and  {d(i)j /dn] 
is  the  vector  of  incident  wave  normal  velocities.  The  above  equation  is  based  on 
sources  distributed  over  the  entire  submerged  portion  of  the  hull  (i.e.,  both  port  and 
starboard  sides).  The  incident  wave  potential  is  given  by: 


i  g  a 

Ui 


exp  [—i  kj  (a:  cos  f3  —  ysin/3)]  exp (kj  zwi ) 


(64) 


Equation  (64)  is  based  on  the  convention  of  the  wave  crest  being  at  the  x  —  y  origin 
at  time  t  =  0.  Fluid  velocities  are  given  by  the  following  derivatives: 


d(pr 

dx 


a  ujj  cos/3  exp  \—ikj  ( xcos/3  —  ysin/3)]  exp (kj  zwi)  (65) 
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d(pr 

dy 

d(f>i 

dz 


— aui  sin  (3  exp  [— i  kj  (a;  cos  (3  —  ?/sin/3)]  exp (kj  zwi ) 
iaui  exp  \—ikj  ( xcos/3  —  ?/sin/3)]  exp  (kj  zwi) 


The  incident  wave  elevation  is  given  by: 


(66) 

(67) 


(j  =  a  exp  [— i  kj  ( xcos/3  —  ?/sin/3)] 


(68) 


When  evaluating  diffraction  potentials,  it  was  found  that  a  significant  amount  of 
time  was  often  required  for  solving  diffraction  source  strengths  in  Equation  (63). 
The  required  solution  time  was  noticeably  greater  than  for  the  solution  of  radia¬ 
tion  source  strengths,  which  utilized  hull  symmetry.  Both  radiation  and  diffraction 
source  strengths  are  solved  using  components  of  the  linear  algebra  library  LA- 
PACK  [15].  The  computational  time  for  solving  a  system  of  equations  with  LA- 
PACK  is  approximately  proportional  to  N3,  where  N  is  the  number  of  unknowns. 
To  improve  computational  speed  for  evaluation  of  diffraction  potentials,  the  inci¬ 
dent  wave  potentials  and  corresponding  diffraction  potentials  can  be  divided  into 
components  which  are  even  and  odd  functions  of  y,  thus  taking  advantage  of  ship 
symmetry.  The  even  and  odd  components  of  the  incident  wave  potential  and  their 
derivatives  are: 


— —  exp  (—ikjx  cos/5)  exp  {ki  zwi)  cos  (kj  y  sin/3)  (69) 

Ui 

aui  cos/3  exp(— %  kj  x  cos/3)  exp  {ki  zwi)  cos  {ki  y  sin/3)  (70) 
—aui  sin/3  exp(— i  kj  x  cos/3)  exp  (kj  zwi )  cos(kj  y  sin/3)  (71) 

i  aui  exp  (— i  ki  xcos/3 )  exp  {ki  zu,i)  cos(A;/  y  sin/3)  (72) 

— —  exp(— /  ki  x  cos/3)  exp  (kj  zwi)  isin(ki  y  sin/3)  (73) 

u>i 

a  ui  cos/3  exp(— /  kj  x  cos/3)  exp  {ki  zwi )  is\n(ki  y  sin/3)  (74) 
— a  ui  sin/3  exp(— /  kj  x  cos/3)  exp  {ki  zwi )  isin(ki  y  sin/3)(75) 
i  a  ui  exp  (— d  ki  x  cos  / 3 )  i  sin  [ki  y  sin  /3)  (76) 


Using  the  even  and  odd  components  of  incident  wave  potential,  source  strengths 
only  have  to  be  solved  for  the  port  side  of  the  hull  using  the  following: 


[. Deven }  {oe£en}  =  {~dtfven/dn}  (77) 

[Dodd]  {a°Jd}  =  {-d^/dn}  (78) 
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where  {a^' en}  is  the  vector  of  source  strengths  on  the  port  side  for  the  even  compo¬ 
nent  of  diffraction  potential,  {d(j)f’en/dn}  is  the  vector  of  incident  normal  velocity 
on  the  port  side  due  to  the  even  component  of  the  incident  wave  potential,  {  a'pd  }  is 
the  vector  of  source  strengths  on  the  port  side  for  the  odd  component  of  diffraction 
potential,  and  {d<J)°Idd / dn}  is  the  vector  of  incident  normal  velocity  on  the  port  side 
due  to  the  odd  component  of  the  incident  wave  potential.  Once  the  source  strengths 
are  solved,  the  total  diffracted  potential  on  the  hull  can  be  determined  by: 


where  {4>p^rt}  is  the  vector  of  diffraction  potentials  on  the  port  side  of  the  hull, 
{0Dar}  ’s  the  vector  of  diffraction  potentials  on  the  starboard  side  of  the  hull, 
[EPort'P°rt]  is  the  influence  matrix  of  velocity  potentials  on  the  port  side  due  to 
sources  on  the  port  side,  and  [ Eport’star ]  is  the  influence  matrix  of  velocity  poten¬ 
tials  on  the  port  side  due  to  sources  on  the  starboard  side. 

The  new  procedure  for  evaluating  diffraction  source  strengths  requires  approxi¬ 
mately  one  quarter  of  the  computation  time  as  the  original  procedure.  Although 
2  sets  of  source  strengths  are  now  solved  (even  and  odd),  each  set  only  has  to  be 
solved  for  half  of  the  hull. 

5.10  Interpolation  of  Inverse  of  Normal  Velocity 
Influence  Matrix  when  Evaluating  Diffraction 
Potentials 

The  database  of  radiation  and  diffraction  potentials  previously  mentioned  allows 
rapid  evaluation  using  interpolation  of  stored  values.  A  comprehensive  database 
encompassing  full  ranges  of  speeds,  headings,  and  wave  frequencies  can  require 
diffraction  computations  for  several  thousand  sets  of  conditions.  The  previous  sub¬ 
section  noted  that  diffraction  potentials  are  evaluated  using  Equations  (57),  (58) 
and  (79).  When  generating  the  database,  the  number  of  required  Green  function 
evaluations  is  greatly  reduced  using  the  following  interpolation  scheme  presented 
in  Reference  1 : 

[DPort*ort(tUe)]  =  Wi  {u  e^\  +  Wl+1  [D*0**0*  (ue,i+1 )]  (80) 

[DPort, star^ j  =  w.  [Dport,star  _.)  ]  +  ^  [jypart^r^j  ]  (gl) 

[E'Por^ort^  j  =  w.  [EP0rt^rt(uJe,)]  +  Wi+l  (ue,i+1)\  (82) 

[E^star{ue)\  =  Wi  [E^star{uJ\  +  wl+l  [E?ort’star Km)j  (83) 
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where  wl  and  wi+1  are  weights  for  linear  interpolation  and  uej  <  u>e  <  weA+ , . 
When  using  the  above  interpolation  scheme,  a  large  amount  of  computational  time 
can  still  be  required  because  the  diffraction  source  strengths  must  still  be  solved 
(solution  of  a^en  and  that  satisfy  Equations  (57)  and  (58))  for  each  combina¬ 
tion  of  ship  speed,  heading,  and  wave  frequency.  This  solution  time  can  be  greatly 
reduced  by  interpolation  of  the  inverse  of  |  f)even]  and  [Dodd]  as  follows: 

[Deven (Ue)]-1  =  Wi  +  wi+1  [Deven(ueti+1)]~l  (84) 

[D°dd{uje)Yl  =  Wi  [Dodd(ujeti)]~1  +  wi+ 1  [D^cu^Y1  (85) 

Using  the  interpolated  inverses,  the  source  strengths  are  then  quickly  evaluated  us¬ 
ing: 

{a%m}  =  [D^^Y1  {~d(j)Ten/dn}  (86) 

{aodd}  =  [D^YYY1  {-d^Y/dn}  (87) 

The  potential  benefits  of  the  above  interpolation  scheme  become  obvious  when  one 
considers  that  a  representative  database  of  diffraction  forces  could  include  7  ship 
speeds  (0,  5, . . .,  30  knots),  13  sea  directions  (0,  15,  . . .,  180  degrees),  and  37  wave 
frequencies  (0.20,  0.25,  . . .,  2.00  rad/s),  giving  a  total  of  3367  seaway  conditions. 
The  same  database  could  have  60  encounter  frequencies  (0.1,  0.2,  . . .,  6.0  rad/s). 

It  should  be  noted  that  there  is  loss  of  accuracy  when  interpolating  influence  ma¬ 
trices  using  Equations  (80)  to  (83).  The  interpolation  of  inverse  matrices  using 
Equations  (84)  and  (85)  likely  leads  to  additional  loss  of  accuracy.  Losses  in  ac¬ 
curacy  will  likely  be  greatest  in  the  vicinity  of  irregular  frequencies  of  encounter, 
at  which  solved  potentials  can  be  sensitive  to  small  variations  in  influence  matrix 
terms. 


5.1 1  New  Procedure  for  Computation  of  Diffraction 
Forces  from  a  Database 


When  predicting  wave  diffraction  forces  using  previously  computed  values,  the  pro¬ 
cedure  presented  in  Reference  1  uses  interpolation  to  obtain  wave  diffraction  force 
in  unit  amplitude  waves  as  a  function  of  ship  speed,  heading,  and  wave  frequency. 
The  main  disadvantage  of  this  approach  is  that  it  does  not  adequately  model  the 
strong  influence  of  wave  encounter  frequency  on  wave  diffraction  force.  To  obtain 
more  accurate  results,  wave  diffraction  forces  are  now  evaluated  as  follows: 


Up  F 


D,u>e 


+ 


uf; 


D,U ' 


(88) 
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where  F^';UJe  is  the  diffraction  force  component  that  is  proportional  to  wave  en¬ 
counter  frequency  and  F^)U  is  the  diffraction  force  component  that  is  proportional 
to  ship  speed.  Reference  1  gives  the  following  equation  for  diffraction  forces: 


if  = 


P  ^ e 


rSb 


U 

UpD  ~ - o — 

UJe  OX 


rij  dS 


(89) 


where  (pc  is  the  potential  of  the  diffracted  waves.  The  corresponding  components 
for  Equation  (88)  are: 


FD,aJe  =  p  [  i  <pD  nj  dS 

(90) 

Jsb 

T 7D,U  f  d<pD 

0  =  P  /  gx  ",  dS 

(91) 

When  creating  a  database  of  hull  radiation  and  diffraction  forces,  each  of  the  above 
components  in  unit  amplitude  waves  are  stored  for  combinations  of  specified  ship 
speeds,  headings,  and  wave  frequencies.  When  accessing  values  from  the  database, 
the  above  components  are  interpolated  as  functions  of  ship  speed,  heading,  and 
wave  frequency. 


5.12  Faster  Versions  of  LAPACK  Routines  for 
Solution  of  Systems  of  Linear  Equations 

When  solving  radiation  and  diffraction  potentials,  a  large  amount  of  CPU  time 
is  devoted  to  solving  systems  of  linear  equations.  Routines  from  the  LAPACK 
library  [15]  are  used  for  solving  systems  linear  equations.  Previously,  functions 
solve  Jinear_equations  and  inverse  from  the  Python  Numeric  library  [16]  were  used 
for  solving  linear  equations  and  inverting  matrices.  The  Python  package  SciPy, 
available  at  no  cost  from  http://www.scipy.org,  has  versions  of  LAPACK  routines 
optimized  for  PC  hardware.  The  SciPy  routines  solve  and  inv  are  now  used  instead 
of  the  Numeric  routines  solve  Jinear. equations  and  inverse. 

Table  3  gives  representative  CPU  times  for  solving  hull  added  mass  and  damp¬ 
ing  coefficients  in  the  frequency  domain  for  HALILAX.  The  values  in  Table  3 
are  based  on  average  CPU  times  over  a  range  of  encounter  frequencies.  As  ex¬ 
pected,  the  amount  of  CPU  time  for  solution  of  Green  functions  is  approximately 
proportional  to  (Nj’ort)2.  The  CPU  time  required  for  matrix  inversion  or  solution 
of  equations  without  matrix  inversion  is  approximately  proportional  to  ( NP°rt )3, 
which  is  also  expected.  The  optimized  inversion  and  solution  routines  from  SciPy 
are  approximately  three  times  faster  than  the  corresponding  routines  from  Numeric. 
Solution  of  velocities  using  matrix  inversion  takes  approximately  three  times  longer 
than  solution  without  matrix  inversion;  however,  matrix  inversion  can  be  useful  for 
determining  matrix  condition  numbers,  as  discussed  in  Section  5.7. 
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Table  3:  Typical  CPU  Times  (seconds)  Per  Encounter  Frequency  for  Prediction  of 
Hull  Added  Mass  and  Damping  Coefficients  in  Frequency  Domain 

Coarse  mesh 

NP°rt  =  224 

Medium  mesh 

Nport  =  433 

Fine  mesh 
NP°rt  _  1040 

Evaluation  of  Green  functions 
and  derivatives 

30 

110 

637 

Inversion  of  influence  matrices 
for  even  and  odd  modes  using 
Numeric  routine  inverse 

2.1 

17 

250 

Inversion  of  influence  matrices 
for  even  and  odd  modes  using 
SciPy  routine  inv 

0.47 

2.9 

45 

Solution  of  velocity  potentials 
for  even  and  odd  modes  using 
Numeric  routine 
solve -linear  .equations 

0.46 

4.8 

81 

Solution  of  velocity  potentials 
for  even  and  odd  modes  using 
SciPy  routine  solve 

0.16 

1.0 

15 

30 
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6  Hydrodynamic  Forces  in  the  Frequency 
Domain  for  HALIFAX 


To  verify  the  improved  methods  presented  in  the  previous  section,  hydrodynamic 
forces  in  the  frequency  domain  have  been  computed  for  HALIFAX.  The  coarse, 
medium,  and  fine  meshes  presented  in  Section  4  have  been  used.  Table  4  gives  pa¬ 
rameters  for  generated  databases  of  ship  hydrodynamic  coefficients  in  the  frequency 
domain.  To  suppress  irregular  frequency  effects,  limiting  envelopes  on  matrix  con¬ 
dition  numbers  were  provided  as  input.  These  envelopes  are  shown  in  Figures  13 
and  14,  and  were  developed  by  examining  initial  predictions  of  hydrodynamic  co¬ 
efficients  and  matrix  condition  numbers  presented  in  Figure  7  to  12 


Table  4:  Parameters  for  Databases  of  Frequency  Domain  Radiation  and  Diffraction 
Forces 


Wave  encounter  frequencies  u>e  0.1,  0.2,  0.3,  . . .,  6.0  rad/s 
Ship  speeds  U  0,  5,  10,  15,  20,  25,  30  knots 

Relative  sea  directions  /3  30,  150,  165,  180  degrees 

Wave  frequencies  ujj  0.20,  0.25,  0.30,  . . .,  2.0  rad/s 


Table  5  gives  required  computation  times  for  generating  the  databases  of  hydrody¬ 
namic  forces  in  the  frequency  domain.  As  discussed  in  Section  5,  wave  diffraction 
computations  can  be  performed  more  quickly  but  somewhat  less  accurately  using 
interpolation  of  influence  matrix  inverses  according  to  Equations  (84)  and  (85).  The 
computational  times  given  in  Table  5  are  based  on  the  more  accurate  approach  of 
interpolation  of  influence  matrices  rather  than  their  inverses. 


Table  5:  Required  CPU  Time  for  Generating  Database  of  Radiation  and  Diffraction 
Computations 


Mesh 

Coarse,  NP°rt  =  224 
Medium,  NP°rt  =  433 
Fine,  NP°rt  =  1040 


CPU  time 
0.7  hours 
2.6  hours 
16.8  hours 
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Figure  13:  HALIFAX  Limiting  Condition  Numbers  for  Longitudinal  Modes  Surge, 
Heave,  and  Pitch 
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Figure  14:  HALIFAX  Limiting  Condition  Numbers  for  Lateral  Modes  Sway,  Roll,  and 
Yaw 
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6.1  Damping  and  x  Derivative  of  Added  Mass 

Figures  15  to  20  show  damping  and  the  x  derivative  of  added  mass  for  HALIFAX. 
These  two  quantities  are  shown  because  they  are  both  used  for  determining  retar¬ 
dation  functions  in  the  time  domain.  The  plots  of  damping  and  x  derivative  of 
added  mass  have  been  non-dimensionalized  to  indicate  their  relative  magnitudes 
for  a  ship  speed  of  20  knots.  This  speed  was  selected  because  it  is  approaching 
the  maximum  speed  at  which  HALIFAX  would  typically  transit  in  a  moderate  to 
heavy  seaway.  Damping  coefficients  from  the  3  different  meshes  have  very  similar 
values.  The  x  derivatives  of  added  mass  are  more  sensitive  to  mesh  size  because 
they  are  small  quantities  arising  from  fore-aft  asymmetry  of  the  ship  hull.  The  good 
agreement  between  the  medium  and  fine  meshes  indicates  numerical  convergence 
with  increasing  number  of  panels. 
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Figure  15:  Dimensionless  Forces  from  Zero  Speed  Surge  Damping  and  x 
Derivative  of  Added  Mass  for  HALIFAX  at  20  Knots 
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Figure  16:  Dimensionless  Forces  from  Zero  Speed  Sway  Damping  and  x 
Derivative  of  Added  Mass  for  HALIFAX  at  20  Knots 
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Figure  17:  Dimensionless  Forces  from  Zero  Speed  Heave  Damping  and  x 
Derivative  of  Added  Mass  for  HALIFAX  at  20  Knots 
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Figure  18:  Dimensionless  Forces  from  Zero  Speed  Roll  Damping  and  x  Derivative 
of  Added  Mass  for  HALIFAX  at  20  Knots 
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Figure  19:  Dimensionless  Forces  from  Zero  Speed  Pitch  Damping  and  x 
Derivative  of  Added  Mass  for  HALIFAX  at  20  Knots 
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Figure  20:  Dimensionless  Forces  from  Zero  Speed  Yaw  Damping  and  x  Derivative 
of  Added  Mass  for  HALIFAX  at  20  Knots 
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6.2  Wave  Excitation  Forces  on  HALIFAX 

Figures  21  and  22  show  wave  excitation  forces  on  HALIFAX  for  a  ship  speed  of  20 
knots  in  stern  quartering  and  bow  quartering  seas.  All  results  are  for  the  medium 
mesh  using  the  following  three  computational  methods,  listing  in  order  of  decreas¬ 
ing  computational  time: 

•  normal  velocity  influence  matrix  [D(tue)\  computed  directly  at  frequency  c oe, 

•  normal  velocity  influence  matrix  [D(ue)\  obtained  using  interpolation  from  ad¬ 
jacent  encounter  frequencies  in  database  (Equations  (80)  and  (81)), 

•  inverse  of  normal  velocity  influence  matrix  [D(ue)\ obtained  using  interpola¬ 
tion  from  adjacent  encounter  frequencies  in  database  (Equations  (84)  and  (85)). 

For  the  most  part,  the  three  approaches  give  excellent  agreement.  Discrepancies 
are  limited  to  when  the  matrix  [D{u>e)\ is  interpolated  for  bow  quartering  seas 
with  wave  frequencies  in  the  vicinity  of  1.0- 1.2  rad/s.  The  associated  encounter 
frequencies  are  of  the  order  2-3  rad/s,  where  the  first  irregular  frequencies  occur 
for  longitudinal  and  lateral  modes,  and  small  errors  in  influence  matrix  terms  can 
lead  to  large  errors  in  computed  solutions.  In  summary,  the  above  results  indicate 
that  reliable  diffraction  forces  can  be  obtained  when  interpolating  the  influence 
matrix  [D]  as  a  function  of  encounter  frequency.  If  the  inverse  matrix  [D]  1  is 
interpolated  as  a  function  of  encounter  frequency,  noticeable  inaccuracies  can  occur 
in  the  vicinity  of  irregular  frequencies. 
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Figure  21:  Wave  Excitation  Forces  for  HALIFAX,  Medium  Mesh,  20  knots,  Stern 
Quartering  Seas  at  30  degrees 
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Figure  22:  Wave  Excitation  Forces  for  HALIFAX,  Medium  Mesh,  20  knots,  Bow 
Quartering  Seas  at  150  degrees 
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6.3  HALIFAX  Motions  in  the  Frequency  Domain 
Using  New  Computational  Methods 

Motions  for  HALIFAX  in  stem  quartering  and  bow  quartering  seas  are  presented 
in  Figures  23  and  24.  The  motions  based  on  directly  computed  hydrodynamic  co¬ 
efficients  agree  very  closely  with  those  from  the  database  with  interpolation  of  in¬ 
fluence  matrix  [D]  as  a  function  of  encounter  frequency.  As  was  observed  for  wave 
excitation  forces,  some  discrepancies  occur  in  the  vicinity  of  irregular  frequencies 
when  the  inverse  matrix  [-D]-1  is  interpolated  as  a  function  of  encounter  frequency. 
When  high  accuracy  is  required  from  a  database  of  diffraction  forces,  the  results 
indicate  that  one  should  interpolate  the  influence  matrix  [D]  rather  than  its  inverse 
[.D]-1  as  a  function  of  encounter  frequency 

For  the  motions  in  stem  quartering  seas  in  Figure  23,  wave  frequencies  have  been 
excluded  that  have  encounter  frequencies  less  than  0.1  rad/s.  At  wave  frequen¬ 
cies  with  encounter  frequencies  below  0.1  rad/s,  predicted  motions  in  unrestrained 
modes  (surge,  sway,  and  yaw)  are  unrealistically  large  due  to  the  absence  of  re¬ 
straining  forces  and  very  small  inertia  and  damping  forces.  For  practical  com¬ 
putations,  this  problem  can  be  avoided  by  a  slight  shift  of  wave  frequency  or  by 
including  viscous  damping  terms. 
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Figure  23:  Frequency  Domain  Motions  for  HALIFAX  from  Direct  and  Database 
Computations,  Medium  Mesh,  20  knots,  Stern  Quartering  Seas  at  30  degrees 
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Figure  24:  Frequency  Domain  Motions  for  HALIFAX  from  Direct  and  Database 
Computations,  Medium  Mesh,  20  knots,  Bow  Quartering  Seas  at  150  degrees 


DRDC  Atlantic  TM  2003-1 04 


45 


7  Hydrodynamic  Coefficients  in  the  Time 
Domain  for  HALIFAX 


Retardation  functions  presented  in  Section  3  have  been  computed  for  HALIFAX 
using  the  frequency  domain  coefficients  presented  in  the  previous  section.  The 
computed  retardation  functions  are  based  on  numerical  integration  of  frequency 
domain  coefficients  for  an  encounter  frequency  range  of  0  to  6  rad/s  with  an  incre¬ 
ment  of  0.05  rad/s.  At  encounter  frequencies  beyond  6  rad/s,  damping  and  the  x 
derivative  of  added  mass  are  assumed  to  be  equal  to  their  infinite  frequency  val¬ 
ues.  Figures  25  to  30  show  the  resulting  speed-independent  and  speed-dependent 
terms  of  retardation  functions.  As  was  done  for  frequency  domain  coefficients  in 
Section  6,  speed-dependent  terms  are  given  for  a  ship  speed  of  20  knots  to  indicate 
their  magnitude  relative  to  speed-independent  terms. 

Figures  25  to  30  show  that  the  speed-independent  terms  are  much  greater  than  the 
speed-dependent  terms  for  HALIFAX  travelling  at  20  knots.  Similar  behaviour 
was  observed  for  frequency  domain  coefficients  in  Section  6.  The  medium  and  fine 
meshes  give  similar  results;  however,  the  coarse  mesh  gives  significant  differences 
in  speed-dependent  terms.  These  differences  can  be  attributed  to  the  relatively  small 
magnitude  of  speed-dependent  terms  and  their  sensitivity  to  computed  x  derivatives 
of  velocity  potentials. 
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Figure  25:  HALIFAX  Surge  Retardation  Contributions  K\\  andUK[[  for  U  =  20 
knots 
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Figure  26:  HALIFAX  Sway  Retardation  Contributions  K22  and  U K22  for  U  =  20 
knots 
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Figure  27:  HALIFAX  Heave  Retardation  Contributions  K-L  and  U K’l,  for  U  =  20 
knots 
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Figure  28:  HALIFAX  Roll  Retardation  Contributions  I<%  and  U K\ \  for  U  =  20 
knots 
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Figure  29:  HALIFAX  Pitch  Retardation  Contributions  K®5,  UK^5,  andU2K §5U 
U  =20  knots 
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Figure  30:  HALIFAX  Yaw  Retardation  Contributions  Kq6,  U K%,  and  U2 K™  for  U 
=  20  knots 
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8  Time  Domain  Simulation  of  Ship  Motions 


The  hydrodynamic  forces  presented  in  this  report  can  form  the  basis  for  simulat¬ 
ing  motions  of  an  unappended  ship.  Ship  motions  can  be  computed  using  a  time 
stepping  procedure  using  the  following  derived  from  Equation  (30): 


([M]  +  [A(U,  oo)])  (r/(f)} 


-  [6(t/)l  {»)(«)} 

-  [  [K(U,t-  dr 

J  — OO 

-  ([Cl  +  [c(u) ])  Mt)} 

+  (F7(t)  +  FD(t)}  (92) 


The  above  system  of  equations  enables  computations  of  accelerations  {?)(£)}  pro¬ 
vided  that  displacements  and  velocities  are  know  up  to  time  t.  Hombeck  [17]  dis¬ 
cusses  several  numerical  methods  that  can  be  used  for  integration  of  ship  motions  in 
the  time  domain.  When  selecting  a  suitable  method,  it  is  desirable  to  minimize  the 
number  of  time  steps  at  which  forces  need  to  be  evaluated  because  of  computational 
requirements.  Another  consideration  is  that  when  integrating  displacements,  both 
first  and  second  derivatives  (velocity  and  acceleration)  are  available.  It  was  decided 
to  use  Taylor  series  for  time-stepping  of  displacements  and  velocities.  The  Taylor 
series  equations  include  contributions  due  to  the  time  derivatives  of  accelerations, 
which  are  estimated  using  a  backward  difference  approach.  The  following  equa¬ 
tions  are  evaluated  in  order  to  obtain  displacements  and  velocities  at  time  t  +  At 
based  on  motions  evaluated  up  to  time  t: 


dt 

r)j(t  +  At) 
r)j(t  +  At) 


23  fij(t) 


Vj(t)  + 
Vj(t)  + 


—  16  rjjit  —  At)  +  5  rjjit 

m 


At  rjj{t)  + 
At  rjjit)  + 


(At)2 

2 

(A  t)2 


ijjit)  + 

dyjjt) 

dt 


-  2A t) 

(At)3  dr) jit) 
6  dt 


(93) 

(94) 

(95) 


When  originally  implementing  the  above  equations  into  the  Python  computer  code, 
arbitrary  time  step  sizes  At  were  facilitated  using  linear  interpolation  of  results 
from  previous  time  steps.  Subsequent  profiling  of  code  execution  revealed  that  a 
large  amount  of  CPU  was  being  devoted  to  interpolation  of  the  retardation  func¬ 
tions  in  Equation  (92).  To  improve  code  efficiency,  the  code  was  re-written  to 
restrict  time  step  size  to  a  single  value  input  at  the  beginning  of  a  simulation.  This 
scheme  gives  much  greater  efficiency  because  retardation  functions  only  need  to 
be  interpolated  at  the  very  beginning  of  a  simulation  using  the  input  time  step  size. 
This  coding  change  has  increased  execution  speed  by  a  factor  of  approximately  30 
for  linear  time  domain  computations. 
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9  Comparison  of  Motions  in  the  Time  and 
Frequency  Domains 


To  verify  hydrodynamic  coefficients  and  derived  motions  predicted  in  the  time  do¬ 
main,  comparisons  have  been  made  with  frequency  domain  motion  predictions. 
The  predictions  are  for  the  HALIFAX  using  the  medium  mesh,  with  433  panels  on 
the  port  side  of  the  hull.  To  give  realistic  roll  motions,  supplementary  roll  damp¬ 
ing  of  2  x  107  Nm/(rad/s)  has  been  used,  which  is  representative  of  damping  from 
appendages  and  viscous  effects.  Similarly,  supplementary  yaw  stiffness  of  2  x  109 
Nm/rad  has  been  used,  which  is  representative  of  the  contributions  of  the  skeg  and 
rudder  to  coursekeeping  stability. 

A  ship  speed  of  20  knots  was  selected  for  comparing  frequency  and  time  domain 
predictions.  This  speed  was  considered  sufficiently  high  to  verify  consistency  of 
speed  dependent  terms.  The  time  domain  computations  are  based  on  a  time  step 
size  of  0.1  s,  which  gives  excellent  numerical  stability.  On  a  PC  computer  with  a 
800  MHz  Pentium  III  processor,  the  time  domain  computations  run  approximately 
30  times  faster  than  real  time.  Figures  3 1  and  32  show  excellent  agreement  between 
predictions  in  the  frequency  and  time  domains  for  both  stem  quartering  and  bow 
quartering  seas.  Note  that  the  present  theoretical  approach  deteriorates  at  very  low 
encounter  frequencies;  thus,  the  predictions  are  limited  to  conditions  with  encounter 
frequencies  >0.1  rad/s.  This  limitation  can  be  handled  from  a  practical  standpoint 
by  slightly  shifting  wave  frequencies  when  travelling  with  forward  speed  in  waves 
aft  of  the  beam. 
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Figure  31:  Motions  for  HALIFAX  from  Frequency  Domain  and  Time  Domain 
Computations,  Medium  Mesh,  20  knots,  Stern  Quartering  Seas  at  30  degrees 


DRDC  Atlantic  TM  2003-1 04 


55 


Time  domain 


Frequency  domain 


Surge 

RAO 

M 


Wave  Frequency  uij  (rad/s) 


Wave  Frequency  ujj  (rad/s) 


Figure  32:  Motions  for  HALIFAX  from  Frequency  Domain  and  Time  Domain 
Computations,  Medium  Mesh,  20  knots,  Bow  Quartering  Seas  at  150  degrees 
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10  Extension  to  Nonlinear  Forces  from 
Incident  Waves  and  Buoyancy 


Time  domain  predictions  presented  thus  far  are  based  on  linearized  forces.  It  is 
relatively  simple  to  extend  the  work  to  include  nonlinear  evaluation  of  forces  due  to 
incident  waves  and  buoyancy.  Furthermore,  these  two  force  components  have  large 
force  magnitudes  relative  to  other  components  such  as  radiation  forces. 

For  evaluation  of  nonlinear  forces,  it  is  necessary  to  introduce  a  new  ship-based 
coordinate  system  xs,  ys,  zs  that  rotates  with  the  ship,  as  shown  in  Figure  33  The 
origin  of  this  coordinate  system  is  at  the  ship  centre  of  gravity.  When  wave-induced 
ship  motions  are  zero,  the  ship-based  coordinate  system  is  identical  to  the  coordi¬ 
nate  system  x,  y,  z  based  on  translating  earth  axes  presented  in  Figure  1. 


Figure  33:  Ship  Referenced  Coordinate  System  for  Large  Angular  Motions 

Force  components  in  the  ship-based  coordinate  system  are  easily  evaluated  as  fol¬ 
lows: 


Fj’(t)  —  —  f  pit)  rij  dS  for  j  =  1  —  6  (96) 

Jsb 

where  nt-  is  the  outward  normal  component  for  mode  j  in  the  ship-based  coordi¬ 
nate  system.  The  following  equation  can  be  used  to  obtain  force  components  in 
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translating  earth  axes  from  force  components  computed  using  ship-based  axes: 

Fi(t)  =  F?(t)  cos 7^6  cos  7/5  -  F2\t)  sin  r]6  cos 775  +  F3s(f)  sin  775  (97) 

F2(f)  =  (cos  ?76  cos  775  sin  774  4-  sin  r)e  cos  774) 

+  F2s(f)  (—sin T/e  sin 775  sin 774  +cos/76  cos 7/4) 

—  F39(f)  sin  774  cos  775  (98) 

F3(f)  =  Ff(t)  (—cos '776  sin  775  cos  774  +  sin  776  sin  774) 

+  F,9  (f )  (sin  ?76  sin  775  cos  774  +  cos  sin  -774) 

+  F9  (f )  cos  774  cos  775  (99) 

The  above  equations  are  based  on  rotations  in  the  order  of  roll,  pitch,  and  yaw 
when  going  from  the  translating  earth  to  ship-based  coordinate  systems.  Similarly, 
the  location  of  a  point  on  the  ship  can  be  converted  to  a  location  in  translating  earth 
axes  as  follows: 


x(t )  =  771  +  xs  cos '776  cos '775  —  ys  sin  ?76  cos  775+  zs  sin  775  (100) 

y(t)  —  rh  +  xs  (cos  776  cos  775  sin  774  +  sin  rjG  cos  774) 

+  ys  (—sin  776  sin  775  sin  774  +  cos  7]q  cos  774) 

—  zs  sin  774  cos  775  (101) 

z(t )  =  773  +  a;5  (—cos '776  sin '775  cos '774  +  sin r]6  sin 774) 

+  y  (sin  776  sin  775  cos  '774  +  cos  776  sin  774) 

+  (f)  cos  774  cos '775  (102) 

^w(f)  =  ^(f)  +  z^f  (103) 

Equations  (102)  and  (103)  are  required  for  evaluating  the  location  of  points  on  the 
hull  within  the  wave  field.  Recall  that  is  the  height  of  the  ship  centre  of  gravity 
above  the  water  surface  when  the  ship  is  at  rest  in  calm  water. 


The  incident  wave  elevation  in  the  time  domain  is  given  by  Equation  (46).  The 
wave  potential  and  its  derivatives  in  the  time  domain  are  given  by: 


$7 

<9$/ 

dx 

dy 

d$! 

dz 


- —  sin  [kj  (a;  cos  f3  —  ysm[3)  —  c oet  —  ej\  exp (kj  zwi)  (104) 

LUi 

a  u>i  cos  / 3  cos  [kj  ( x  cos  (3  —  y  sin  (3)  —  a )et  —  e/] 
x  exp (&/  zwi)  (105) 

—a  ujj  sin  (3  cos  [ kj  ( x  cos  (3  —  y  sin  (3)  —  c oet  —  ej] 
x  exp  (A;/  zwi)  (106) 

a  ujj  sin  [kj  (a;  cos  (3  —  ysin/3)  —  c oet  —  ej\  exp  (kj  zwi)  (107) 
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~dt 


(108) 


=  —  - — c oe  cos  [ki  (x  cos (3  —  y  sin/3) 

Ui 

x  exp (k!  zwi) 


ujf,  t 


ei] 


Using  the  above  terms,  the  combined  pressure  due  to  incident  wave  forces  and 
buoyancy  can  be  evaluated  as  follows: 


fd§j  d  <f>r  \ 

pit)  =  -p  f  — - U  +  g  zwA  for  zwi  <  (  (109) 

p(t)  =  0  for  zwi  >  C  (110) 


Figures  34  to  36  show  results  of  computations  for  verifying  nonlinear  buoyancy 
force  predictions.  For  small  displacements,  the  nonlinear  buoyancy  predictions 
show  close  agreement  with  linear  buoyancy  forces.  For  large  heave  displacements, 
nonlinear  buoyancy  forces  approach  limiting  values  as  the  ship  hull  becomes  com¬ 
pletely  submerged  (negative  heave)  or  emerged  (positive  heave).  The  nonlinear 
roll  buoyancy  moments  behave  as  expected,  with  symmetry  about  zero  roll  angle 
and  vanishing  righting  moment  at  large  roll  angles.  The  nonlinear  pitch  buoyancy 
moments  behave  in  a  manner  similar  to  the  heave  forces,  with  slope  magnitude 
reducing  when  portions  of  the  hull  become  completely  submerged  or  emerged. 
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Figure  34:  Heave  Buoyancy  Forces  for  HALIFAX,  Medium  Mesh 
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Figure  35:  Roll  Buoyancy  Moments  for  HALIFAX,  Medium  Mesh 
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Figure  36:  Pitch  Buoyancy  Forces  for  HALIFAX,  Medium  Mesh 
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11  Comparison  of  Numerical  Predictions 
and  Experimental  Results  for  Heave  and 
Pitch 


The  software  testing  presented  thus  far  has  been  limited  to  verification,  including 
comparison  of  predictions  in  the  frequency  and  time  domains.  Validation  of  ship 
motion  predictions  generally  requires  comparison  with  measurements  from  model 
tests  and/or  sea  trials.  The  current  motion  predictions  do  not  include  viscous  or 
appendage  forces,  which  can  significantly  influence  lateral  plane  motions;  however, 
predicted  vertical  plane  predictions  can  be  validated  using  experimental  data. 

Figures  37  and  38  show  comparisons  of  heave  and  pitch  predictions  with  exper¬ 
imental  results  for  the  CPF  hydroelastic  model  in  head  seas  [7].  The  numerical 
predictions  include  the  influence  of  restrained  surge  motions.  The  numerical  pre¬ 
dictions  show  very  good  agreement  with  experimental  motions.  As  was  observed 
in  Section  9,  frequency  domain  predictions  correspond  very  closely  to  linear  time 
domain  predictions.  The  nonlinear  predictions,  which  are  for  a  wave  steepness 
H/ A  of  1/20,  vary  noticeably  from  the  linear  predictions.  For  the  current  heave  and 
pitch  cases,  it  appears  that  including  nonlinear  buoyancy  and  incident  wave  forces 
does  not  increase  the  accuracy  of  predictions.  For  intermediate  wave  frequencies 
at  a  Froude  number  of  0.20,  Figure  38  suggests  that  including  nonlinearities  can 
actually  give  less  accurate  heave  predictions.  This  observation  could  be  due  to 
modelling  of  larger  amplitude  waves  using  linear  wave  theory.  To  improve  accu¬ 
racy  of  nonlinear  predictions  in  regular  waves,  a  new  regular  wave  model  will  be 
implemented  using  Stokes  second  order  wave  theory  (see  Sarpkaya  and  Isaacson 
[18]). 

As  mentioned  in  Section  9,  the  linear  time  domain  simulations  run  approximately 
30  times  faster  than  real  time  when  simulating  a  full-scale  ship  using  a  time  step 
size  of  0.1  s.  When  nonlinear  buoyancy  and  incident  wave  forces  are  included  using 
the  medium  mesh,  computational  speed  is  reduced  to  approximately  5  times  faster 
than  real  time. 
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Experiments 
•  H/X  =  1/30 
A  H/X  =  1/20 
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Figure  37:  Heave  and  Pitch  for  CPF  Hydroelastic  Model  in  Head  Seas,  Fn  =  0. 12 
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Experiments 
•  H/X  =  1/30 
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Figure  38:  Heave  and  Pitch  for  CPF  Hydroelastic  Model  in  Head  Seas,  Fn  =  0.20 
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12  Conclusions 


New  software  has  been  developed  for  predicting  hydrodynamic  forces  and  result¬ 
ing  motions  in  the  time  domain  for  an  unappended  ship  hull  in  waves.  The  ship  is 
assumed  to  travel  at  steady  speed  and  heading  relative  to  the  waves.  The  hydrody¬ 
namic  forces  in  the  time  domain  are  evaluated  using  previously  computed  values 
in  the  frequency  domain.  Several  improvements  have  been  made  to  the  accuracy, 
efficiency,  and  robustness  of  hydrodynamic  force  predictions  in  the  frequency  do¬ 
main.  Comparisons  between  predicted  motions  in  the  frequency  and  time  domains 
indicate  consistency  between  the  two  approaches.  Initial  validation  using  experi¬ 
mental  motions  for  heave  and  pitch  indicates  that  predictions  in  the  frequency  and 
time  domains  are  giving  good  results.  Nonlinear  effects  arising  from  incident  wave 
and  buoyancy  forces  can  be  included  in  time  domain  predictions. 

In  the  near  term,  forces  from  rudders,  bilge  keels,  and  other  appendages  will  be 
incorporated  into  numerical  predictions  of  ship  motions  in  waves.  In  the  longer 
term,  the  motions  of  a  freely  maneuvering  ship  in  waves  will  be  predicted. 
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Symbols 
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a 

B 

[B] 

[b] 
[G] 

[c] 

CB 

[D] 


even 


Fn 

FP 


Go 


9 

[H\ 

[K] 

[K°] 

[Ku] 

[Kuu] 

KG 

ki 

L 

LCB 

[M\ 

nk 

nk 

ATport 

P 

odd 

P 

port 


added  mass  matrix 
wave  amplitude 
ship  beam 

frequency  domain  damping  matrix 

dynamic  damping  matrix 

hydrostatic  stiffness  matrix 

dynamic  stiffness  matrix 

centre  of  buoyancy 

normal  flow  velocity  influence  matrix 

superscript  denoting  velocity  potential  is  even  function  of  y 

diffracted  wave  force  vector 

diffracted  force  component  proportional  to  U  for  mode  j 

diffracted  force  component  proportional  to  uje  for  mode  j 

incident  wave  force  vector 

force  for  mode  j  in  ship-based  coordinates 

Froude  number 

fore  perpendicular 

frequency  dependent  term  of  Green  function  relative  to  zero 
frequency  Green  function 

frequency  dependent  term  of  Green  function  relative  to  infinite 
frequency  Green  function 
gravitational  acceleration 
x  velocity  influence  matrix 
retardation  function  matrix 

zero-speed  component  of  retardation  function  matrix 

component  of  retardation  function  matrix  proportional  to  U 

component  of  retardation  function  matrix  proportional  to  U2 

height  of  centre  of  gravity  above  baseline 

incident  wavenumber 

ship  length  between  perpendiculars 

longitudinal  centre  of  buoyancy 

inertial  matrix 

outward  normal  component  k  for  translating  earth  axes 
outward  normal  component  k  for  ship-based  axes 
number  of  panels  on  port  side  of  hull 
superscript  denoting  velocity  potential  is  odd  function  of  y 
pressure 

superscript  denoting  port  side  of  hull 
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distance  from  field  point  to  source 

distance  from  field  point  to  source  image 

roll  radius  of  gyration 

pitch  radius  of  gyration 

yaw  radius  of  gyration 

surface  of  floating  body 

surface  of  hull  panel  k 

superscript  denoting  starboard  side  of  hull 

draft  at  midships 

trim  by  stern 

ship  speed 

interpolation  weight  factor 
translating  earth  axis  coordinates 
ship-based  axis  coordinates 
source  location 

height  of  ship  centre  of  gravity  above  water  surface  in  calm  water 

z  coordinate  relative  to  calm  water  surface 

wave  direction  relative  to  ship 

time  step  size 

Kroenecker  delta  function 

phase  lead  of  incident  wave 

incident  wave  elevation 

motion  displacement  in  mode  j 

matrix  condition  number 

water  density 

source  strength  vector 

delay  time  for  convolution  integral 

incident  wave  potential  in  time  domain 

diffracted  wave  potential  in  frequency  domain 

incident  wave  potential  in  frequency  domain 

radiated  wave  potential  due  to  motion  mode  k 

normalized  potential  for  mode  k  at  zero  or  infinite  frequency  limit 

incident  wave  frequency 

wave  encounter  frequency 

transition  encounter  frequency  for  Green  function  computation 
ship  mass  displacement 
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